Concentration of Antioxidant Compounds from Calendula officinalis through Sustainable Supercritical Technologies, and Computational Study of Their Permeability in Skin for Cosmetic Use

The growing interest in the cosmetic industry in using compounds of natural and sustainable origin that are safe for humans is encouraging the development of processes that can satisfy these needs. Chlorogenic acid (CHA), caffeic acid (CAF) and ferulic acid (FA) are three compounds widely used within the cosmetic industry due to their functionalities as antioxidants, collagen modifiers or even as radiation protectors. In this work, two advanced separation techniques with supercritical CO2 are used to obtain these three compounds from Calendula officinalis, and these are then evaluated using a computational skin permeability model. This model is encompassed by the COSMO-RS model, the calculations of which make it possible to study the behaviour of the compounds in the epidermis. The results show that both CAF and FA are retained in the stratum corneum, while CHA manages to penetrate to the stratum spinosum. These compounds were concentrated by antisolvent fractionation with super-critical CO2 using a Response Surface Methodology to study the effect of pressure and CO2 flow rate. CHA, CAF and FA were completely retained in the precipitation vessel, with concentrations between 40% and 70% greater than in the original extract. The conditions predicted that the optimal overall yield and enrichment achieved would be 153 bar and 42 g/min.


Introduction
In recent years, there has been great interest in the cosmetic industry, among others, in including compounds of natural and sustainable origin in the formulation of products. This is due to, on the one hand, growing general concern related to the possible harmful effects of synthetic compounds and, on the other hand, awareness regarding the responsible consumption of these sustainable products, causing an increase in their demand. Additionally, the traditional use of plants as natural remedies, or even as cosmetics, makes them ideal candidates from which to obtain compounds that meet these requirements.
Calendula officinalis L., also known as marigold, is an aromatic herbaceous plant that belongs to the genus Calendula (Asteraceae), native to the Mediterranean countries. This genus includes approximately 25 species, of which C. officinalis is the only one used extensively clinically worldwide [1], traditionally as a skin remedy for dermatological problems, such as inflamed skin, redness, minor burns, or ulcers, as well as for acne or fungal eruptions [2]. Over the years, studies have demonstrated these properties, which have been duly collected in various monographs such as those by the European Scientific Cooperative On Phytotherapy (ESCOP), the European Medicines Agency (EMA) or the World Health Organization (WHO), where wound healing and anti-inflammatory actions stand Taking into account the above, one of the objectives of this work is to obtain, in a sustainable way and by combining two advanced separation techniques, an extract enriched in antioxidant compounds from Calendula officinalis. To do this, first, an extraction with supercritical CO 2 will be carried out to defat the plant material and facilitate obtaining the antioxidant substances. After maceration in ethanol, a study of the influence of pressure (80-160 bar) and of the CO 2 flow rate (10-60 g/min) on the SAF process will be performed. This will be performed using a Response Surface Methodology (RSM) based on Central Composite Design (CCD), which will provide the values for the studied magnitudes and the sequence to follow when doing the experiments. The compounds chosen to be concentrated on were FA, CHA and CAF, due to their antioxidant properties [21][22][23].
On the other hand, the fact that a compound is of natural or synthetic origin does not affect the possibility that it may have harmful or beneficial effects on people's health, so it is necessary to study its interaction with human tissues. In vivo and in vitro tests are expensive and difficult to perform, in addition to the fact that there are several regulations that strongly recommend or require the use of alternatives to animal studies (EU REACH regulation) or that even prohibit the marketing of cosmetics with ingredients tested on animals (EU Cosmetics Regulation EC 1223/2009) [46]. In this regard, in silico experimentation is of great interest, because it makes the study of the interaction of solutes with biological membranes possible, and therefore, their possible toxicity. There are various predictive models for the bioavailability of solutes based on quantitative structure-permeability relationship (QSPR) models, molecular dynamics (MD) simulation, or mechanistic models derived from first principles such as mass balance, relying on additional assumptions such as Fick's laws of diffusion [47].
In this research framework, in silico research on the interaction of solutes with biomembranes, the study of molecular mechanisms behind skin permeation stands out, since it is a tool that makes it possible to identify inappropriate ingredients that may be dangerous in a cosmetic composition [48,49]. One of the most innovative models is the one proposed by Schwöbel et al. [46] which is based on the use of the Conductor-like Screening model for Real Solvents (COSMO-RS). From the extensions of this model-COSMOperm, COS-MOmic and COSMOplex-a computational skin model that makes it possible to calculate the permeability of different compounds in it is generated. This model was applied to study the interaction of the tracked compounds of C. officinalis, CHA, CAF, FA, with the skin and thus to calculate their permeabilities.

Plant Material
Calendula officinalis L. flowers were collected from an ecological cultivation in Huesca (Spain) by Valentia, an association where people at risk of social exclusion work. The plant material was dried in a dryer at room temperature. Dried Calendula flowers had an 11.29 ± 0.48 wt% humidity content, determined ten times with a Sartorious model MA 40 Moisture Analyzer (Goettingen, Germany). The plant material was ground with an electric grinder and sieved with a vibratory sieve shaker (CISA model BA 300N, Barcelona, Spain). The average diameter of particles was 0.5 mm. This was adjusted to a normal distribution according to ASAE S319.3 from the American National Standards Institute [50].

Supercritical CO 2 Extraction (SCE)
In a first stage, the plant material was defatted with supercritical CO 2 in a Waters laboratory scale plant (SFE-1000F-2-FMC10 System, Pennsylvania, USA). Its scheme has been published previously [51] and its main components consist of an extraction vessel of 1 L (E), and two 0.5 L collectors (C1, C2). All three vessels were jacketed to maintain a constant temperature. The CO 2 stored in a bottle was passed through a cooling bath (CB) to keep it liquid and then propelled through a pump (P2) towards E. A heat exchanger (HE) was used to ensure that the CO 2 passing into E was above the critical temperature. Flow rate and temperature were automatically controlled, as was the pressure in E. The pressure in C1 and C2 was controlled by manual back pressures (MBPR).
The procedure is the same as that performed in previous works [51]. The same proportion of plant material/inert glass bead was used (100 g/200 g), which allows a better CO 2 -solid contact, and therefore facilitates better extraction. The complete extraction process consisted of 4 static-dynamic cycles and the extraction conditions were 350 bar, 40 • C in E, 90 bar, 45 • C in C1 and 30 bar and 30 • C in C2. The CO 2 flow rate was 60 g/min. Once the machine was depressurized, the extracts collected from C1 and C2 and the plant material recovered from E were stored in a freezer for further experiments and analysis.
The supercritical extraction yield (Y SCE ) was calculated through Equation (1): where mass(g) C1 and mass(g) C2 are the mass in grams of the extract deposited in collectors C1 and C2, respectively, of the SCE device, and mass(g) vegetal material is the initial mass in grams of the dried plant material loaded in the extractor.

Maceration and Supercritical Antisolvent Fractionation (SAF) Processes
Maceration with ethanol (EtOH) was performed to obtain the polar and active compounds of C. officinalis flowers. A 3 L volume of the solvent was used to macerate 300 g of vegetal material, previously defatted with SCE, for 48 h at room temperature (25 • C). Then, a rotatory evaporator (Büchi R-200, Flawil, Switzerland) was used to remove the solvent and obtain the dry extract. Equation (2) indicates the formula to calculate the maceration yield: where mass plant extract (g) was the mass of the solvent-free extract obtained and mass vegetal material (g) was the initial mass, in grams, of plant material subjected to maceration and on which the SCE process was applied.
To prepare the feed solution (FS) for each of the SAF experiments, the extract obtained (without solvent) was redissolved in ethanol at 3% (wt%) and filtered through Cellulose Acetate filters with a pore size of 0.22 µm. SAF experiments were carried out by a laboratoryscale plant (Waters, Pennsylvania, USA), the scheme of which is represented in Figure 1. The plant was equipped with a CO 2 pump (P-SCF), an FS pump (P-LIQ), a 0.5 L precipitation vessel (PV), and a 0.5 L downstream vessel (DV). Temperatures and flow rates of both CO 2 and FS were automatically controlled, as was pressure in PV (ABPR), but pressure in DV was controlled by manual back pressure regulator (MBPR) [52].
A typical SAF experiment starts when CO 2 is pumped by P-SCF through a heat exchanger (HE) to PV and the rest of the plant. Once the selected condition of pressure, temperature and CO 2 flow rate in the system have stabilised, FS is pumped into the PV through an injector (nozzle Ø = 100 mm). In this vessel, the insoluble compounds in the supercritical mixture (CO 2 + EtOH) precipitate, while the soluble compounds are gathered in DV. Once the flow of FS has finished, 30 mL of pure ethanol is pumped to draw the remaining FS from the pipes. Finally, a flow of CO 2 is maintained to eliminate the residual ethanol from the solid precipitated in PV. The conditions used were selected on the basis of the previous knowledge of the research group: to avoid thermal degradation, temperature in PV was fixed at 40 • C; to guarantee the supercritical state of the mixture (CO 2 + EtOH), the FS flow rate and the FS concentration were fixed at 0.45 mL/min and 3% (wt%), respectively, for all experiments [53]. Then the pressure in PV and the CO 2 flow rate were studied and varied between 80 and 160 bar and 10 and 60 g/min, respectively. in Figure 1. The plant was equipped with a CO2 pump (P-SCF), an FS pump (P-LIQ), a 0.5 L precipitation vessel (PV), and a 0.5 L downstream vessel (DV). Temperatures and flow rates of both CO2 and FS were automatically controlled, as was pressure in PV (ABPR), but pressure in DV was controlled by manual back pressure regulator (MBPR) [52]. A typical SAF experiment starts when CO2 is pumped by P-SCF through a heat exchanger (HE) to PV and the rest of the plant. Once the selected condition of pressure, temperature and CO2 flow rate in the system have stabilised, FS is pumped into the PV through an injector (nozzle Ø = 100 mm). In this vessel, the insoluble compounds in the supercritical mixture (CO2 + EtOH) precipitate, while the soluble compounds are gathered in DV. Once the flow of FS has finished, 30 mL of pure ethanol is pumped to draw the remaining FS from the pipes. Finally, a flow of CO2 is maintained to eliminate the residual ethanol from the solid precipitated in PV. The conditions used were selected on the basis of the previous knowledge of the research group: to avoid thermal degradation, temperature in PV was fixed at 40 °C; to guarantee the supercritical state of the mixture (CO2 + EtOH), the FS flow rate and the FS concentration were fixed at 0.45 mL/min and 3% (wt%), respectively, for all experiments [53]. Then the pressure in PV and the CO2 flow rate were studied and varied between 80 and 160 bar and 10 and 60 g/min, respectively.
Equation (3) was used to calculate mass recovery yields for the precipitation (PV) and the downstream (DV) vessel fractions, YPV and YDV, respectively, of the SAF device: (wt%) = mass fraction collected mass of extract in FS · 100 where i refers to the vessel from which the mass is collected: PV or DV. The overall recovery yield of the process, YSAF, was calculated using Equation (4): Equation (3) was used to calculate mass recovery yields for the precipitation (PV) and the downstream (DV) vessel fractions, Y PV and Y DV , respectively, of the SAF device: where i refers to the vessel from which the mass is collected: PV or DV. The overall recovery yield of the process, Y SAF , was calculated using Equation (4):
where i is the considered antioxidant (CHA, CAF, FA) and j refers to the fraction in which the compound was analysed, i.e., FS, PV or DV. Once Ci/j had been obtained, the enrichment ratio / was obtained for each compound using Equation (6): where, again, i refers to the compound analysed (CHA, CAF, FA) and j is PV or DV. This equation can be modified to calculate how much the PV or DV fractions are enriched in Equation (5) was used to calculate the concentration of a compound in each sample analysed: where i is the considered antioxidant (CHA, CAF, FA) and j refers to the fraction in which the compound was analysed, i.e., FS, PV or DV. Once C i/j had been obtained, the enrichment ratio E i/j was obtained for each compound using Equation (6): where, again, i refers to the compound analysed (CHA, CAF, FA) and j is PV or DV. This equation can be modified to calculate how much the PV or DV fractions are enriched in antioxidant compounds. Equation (7) shows how the enrichment ratio of all compounds analysed in PV fraction can be calculated:

Experimental Design and Statistical Analysis
One of the aims in this work is to evaluate the influence of pressure and CO 2 flow rate on the SAF process and then to optimize the conditions to obtain the best results on yield and enrichment ratios. To statistically study this, a response surface methodology (RMS) based on central composite design (CCD) was used. The values were set between 80 and 160 for pressure in PV, X P , and 10-60 g/min for CO 2 flow rate, XQ CO 2 . The rest of the parameters in the experiments-temperature in PV and DV, pressure in DV, and FS flow rate-were constant and were set according to previous experience to ensure the supercritical conditions of the CO 2 -ethanol mixture [51].
The mathematical model of a two variable CCD makes it possible to correlate, through Equation (8), a dependent variable, Y, with some independent variables, X i and X j : where β 0 is a constant coefficient, β i is the linear coefficient, β ii is the quadratic coefficient, and β ij is an interaction coefficient. In this work, the dependant variable refers to both the yields of the SAF process (Y PV , Y DV and Y SAF ) and the enrichment ratios of the tracked compounds (CHA, CAF and FA), while the independent variables, which are the variables under study, are pressure (X P ) and CO 2 flow rate (XQ CO 2 ). The software used to perform the CCD design was Minitab ® 18. This design provides 13 random experiments, 5 of which are replicates of the central conditions according to the range levels of the variables used, shown in Table 1. This software was also used to determinate the values of each coefficient, β, in Equation (8) and their significance (when p < 0.05). Finally, Minitab ® 18 was also used to obtain the optimal conditions for the maximum overall recovery yield (Y SAF ) and maximum enrichment ratio of the bioactive compounds (CHA, CAF, FA) (Y ALL/PV ).

Application of the Skin Model
COSMO-RS is a continuum solvation model that makes it possible to calculate thermodynamic properties, generating a three-dimensional distribution of surface polarization charge-densities, σ, from optimized 3D structures of molecules [53,54]. This σ is used, by means of an efficient statistical thermodynamic model of pairwise molecular surface interactions, to calculate the σ-potential, which gives the chemical potential of a surface segment of polarity in a particular solvent. Then, those segments are added and corrected by a combinatorial term to calculate the chemical potential of a compound in a pure or mixed solvent [53]. This can be applied, for example, to calculate the partition coefficient, K, of a solute between two liquid phases [55].
COSMOperm is an extension of the COSMO-RS model that applies these calculations in inhomogeneous systems such as biomembranes. The COSMOmic method, included within COSMOperm, allows the calculation of partition coefficients of membranes or micelles, chemical potentials, and free energies of solutes in a layered system, which in turn makes it possible to obtain their distribution within the membrane [56][57][58][59][60]. The distribution of the different compounds that make up a biomembrane can be obtained using classic molecular dynamics simulations or using another extension of the COSMO-RS model, called COSMOplex. This extension generates divided liquid layers of variable composition as a representation of a biomembrane from information about its structure [61]. This makes it possible to generate a skin model with which COSMOperm works to calculate the resistance of the membranes, the permeability of individual compounds as well as their position in them, and the permeation pathway [46]. The computational skin model used in this work is based on the model proposed by Schwöbel et al. [46]. This model is based on the division of the outermost layer of the skin, the epidermis, into different compartments, or layers, with their own structure and cellular compositions [47], as can be seen in Table 2. The skin penetration model is based on a series of resistances, R, from which the permeability coefficient, K p , can be calculated [62], as indicated in Equation (9).
where R skin is the overall skin resistance. To calculate R skin , it is first necessary to calculate the sum of all the resistances of each of the compartments that make up it, thus obtaining R stratified-cells as indicated by Equation (10): where R SC , R SG , R SS , and R SB are the resistances of the SC, SG, SS, and SB compartments, described in Table 2, respectively. Starting from R stratified-cells and also taking into account the shunt pathway, the R skin is calculated following Equation (11): where R shunt is the resistance of the shunt pathway and is kept constant 1/R shunt = 2 × 10 −11 m/s. In addition, the transport of the compounds through the compartments is calculated by Equation (12): where i refers to the compartment (SC, SG, SS or SB), R i,trans is the mechanism of transcellular absorption (through keratin-corneocytes by partitioning into and out of the cell membrane), and R i,inter is the mechanism of intercellular absorption (trough corneocytes in the lipid-rich extracellular regions). The skin model was applied by COSMO-RS, using the COSMOtherm software package for each of the tracked compounds in this work: CHA, CAF and FA. First, the pre-optimized 3D chemical structures from PubChem database were obtained, and then those structures were refined by using Gaussian version 9.0 with a DFT parametrization bvp86-TZVP (cartesian coordinates for optimized geometries can be found in the Supplementary Material). The 3D structures of the compounds and their surface charge density can be observed in Figure 3. The skin model, included in COSMOplex module, was used together with the COSMOperm extension to calculate the permeabilities of each compound as well as their position in the skin model and the resistances of each compartment. optimized 3D chemical structures from PubChem database were obtained, and then those structures were refined by using Gaussian version 9.0 with a DFT parametrization bvp86-TZVP (cartesian coordinates for optimized geometries can be found in the Supplementary Material). The 3D structures of the compounds and their surface charge density can be observed in Figure 3. The skin model, included in COSMOplex module, was used together with the COSMOperm extension to calculate the permeabilities of each compound as well as their position in the skin model and the resistances of each compartment.

SCE and Feed Solution Preparation
The pre-treated material, previously ground and sieved, was defatted in a first stage to facilitate the obtaining of the antioxidant fraction. This stage was usually carried out by means of maceration in hexane [41,45], but in a more recent study, supercritical extraction with CO2 was used [51]. This is due to the advantages offered by this type of extraction; since it is a non-toxic solvent, its polarity facilitates the extraction of lipophilic compounds and generates solvent-free final products because it can be easily removed by lowering the pressure. The supercritical extraction yield, YSCE, was 8.3%, which is similar to the yields obtained by D. Baumman et al. The polar fraction of Calendula officinalis was obtained by maceration in ethanol (EtOH) of the plant material after the SCE. The extraction yield obtained with this ethanolic maceration, YEtOH, was 7.2%. Other authors [65] have obtained similar extraction yield (8.0%), although the plant material was not previously defatted and the proportional mass plant (g):EtOH (mL) was 1:6 instead of 1:10, as used in this work.

SAF Yields Statistical Analysis
Experimental values for YPV (wt%), YDV (wt%) and YSAF (wt%) for every experiment are gathered in Table 3. YPV varied from 12.3% (run 10; 120 bar-60 g/min) to 47.0% (run 1; 80 bar-35 g/min), whereas YDV oscillated between 1.3% (run 1; 80 bar-35 g/min) and 42.8% (run 12; 148 bar-53 g/min). It can be seen that at low pressures, PV yields are higher than DV yields, but this changes as the pressure increases; then, the yields equalize until YDV are greater than YPV. The greatest difference was found in run 1, where YPV was 36 times the DV yield. Overall yields, YSAF, varied between 31.1% (run 3; 92 bar-53 g/min) and 73.7% (run 12; 148 bar-53 g/min). As seen previously in other works [51], a full recovery of the entire mass of solutes contained in the feed solution is not possible due to the effect of two

SCE and Feed Solution Preparation
The pre-treated material, previously ground and sieved, was defatted in a first stage to facilitate the obtaining of the antioxidant fraction. This stage was usually carried out by means of maceration in hexane [41,45], but in a more recent study, supercritical extraction with CO 2 was used [51]. This is due to the advantages offered by this type of extraction; since it is a non-toxic solvent, its polarity facilitates the extraction of lipophilic compounds and generates solvent-free final products because it can be easily removed by lowering the pressure. The supercritical extraction yield, Y SCE , was 8.3%, which is similar to the yields obtained by D. Baumman et al. The polar fraction of Calendula officinalis was obtained by maceration in ethanol (EtOH) of the plant material after the SCE. The extraction yield obtained with this ethanolic maceration, Y EtOH , was 7.2%. Other authors [65] have obtained similar extraction yield (8.0%), although the plant material was not previously defatted and the proportional mass plant (g):EtOH (mL) was 1:6 instead of 1:10, as used in this work.

SAF Yields Statistical Analysis
Experimental values for Y PV (wt%), Y DV (wt%) and Y SAF (wt%) for every experiment are gathered in Table 3. Y PV varied from 12.3% (run 10; 120 bar-60 g/min) to 47.0% (run 1; 80 bar-35 g/min), whereas Y DV oscillated between 1.3% (run 1; 80 bar-35 g/min) and 42.8% (run 12; 148 bar-53 g/min). It can be seen that at low pressures, PV yields are higher than DV yields, but this changes as the pressure increases; then, the yields equalize until Y DV are greater than Y PV . The greatest difference was found in run 1, where Y PV was 36 times the DV yield. Overall yields, Y SAF , varied between 31.1% (run 3; 92 bar-53 g/min) and 73.7% (run 12; 148 bar-53 g/min). As seen previously in other works [51], a full recovery of the entire mass of solutes contained in the feed solution is not possible due to the effect of two phenomena: the dragging of the most volatile components through the vent valve [66] and the deposition of materials in dead spaces.
The coefficients of Equation (8) were obtained for Y PV , Y DV and Y SAF and can be found in Table 4, along with their level of significance p and the coefficient R 2 and the deviation s of the fitted mathematical model. The statistical analysis reveals that Y PV depends on all terms (β 1 , β 2 , β 11 and β 12 ) except the quadratic term of CO 2 (β 22 ) and all of them are statistically significant (p < 0.05) except the pressure term (β 1 ). Y DV and Y SAF depend on all terms being statistically significant only for Y SAF . For Y DV , the statistically significant terms are the pressure (β 1 ), the quadratic term of pressure (β 11 ), and the cross term (β 12 ). Table 3. Operational conditions of pressure, X P , and CO 2 flow rate XQ CO 2 , for every run of the CCD design of the SAF process as well as the corresponding results for the yields and enrichment ratios, defined by Equations (3)-(7).

Run
Exp. Run Order X P (bar)   Figure 4, the contour plots corresponding to the surfaces defined by Equation (8) are shown for all yields (Y PV , Y DV and Y SAF as 1a, 1b and 1c, respectively). In Figure 4a, for flow rates above 12 g/min and below 24 g/min, Y PV decreases as pressure increases, while for CO 2 flow rates outside this range, Y PV first decreases and then increases as pressure increases. For a fixed pressure, an increase in CO 2 flow rate causes Y PV to decrease except for high pressures, whereas the CO 2 flow rate increases, Y PV increases. For the studied ranges of pressure and CO 2 flow rates, the highest Y PV is found at low pressures (80-83 bar) and at low CO 2 flow rates (10-16 g/min). Another maximum in Y PV yield could be found at high pressures and CO 2 flow rates (133-160 bar, 51-60 g/min). According to Figure 4b, Y DV increases as pressure increases for a set CO 2 flow rate. On the other hand, for a fixed pressure, a different behaviour is observed if the pressure is lower or higher than 133 bar: for low pressures, as the flow increases, Y DV decreases, whereas the opposite effect is observed at higher pressures. In fact, the highest Y DV is found at high pressures and CO 2 flow rates (133-160 bar, 47-60 g/min). This may be due to a greater solubility of components and a higher dragging effect of the solutes towards DV. Analysing Figure 4c, it can be seen that Y SAF follows a similar behaviour to Y PV , whereby the overall recovery of solutes from the feed solution is minimal at low pressures and high CO 2 flow rates, whereas Y SAF reaches a maximum when the two variables simultaneously reach their maximum (158-160 bar, 59-60 g/min) or minimum (80-83 bar, 10-12 g/min) values in the intervals considered. flow rates (133-160 bar, 47-60 g/min). This may be due to a greater solubility of components and a higher dragging effect of the solutes towards DV. Analysing Figure 4c, it can be seen that YSAF follows a similar behaviour to YPV, whereby the overall recovery of solutes from the feed solution is minimal at low pressures and high CO2 flow rates, whereas YSAF reaches a maximum when the two variables simultaneously reach their maximum (158-160 bar, 59-60 g/min) or minimum (80-83 bar, 10-12 g/min) values in the intervals considered.

Enrichment Ratio Analysis
In addition to the SAF global yields, ferulic acid (FA), caffeic acid (CAF) and chlorogenic acid (CHA) were monitored for the feed solution fraction (FS), precipitation vessel (PV) and downstream vessel (DV). The enrichment ratios Ei/j and EALL/PV, defined by Equations (6) and (7), are gathered in Table 3 for all experiments performed. No Ei/DV is

Enrichment Ratio Analysis
In addition to the SAF global yields, ferulic acid (FA), caffeic acid (CAF) and chlorogenic acid (CHA) were monitored for the feed solution fraction (FS), precipitation vessel (PV) and downstream vessel (DV). The enrichment ratios E i/j and E ALL/PV , defined by Equations (6) and (7), are gathered in Table 3 for all experiments performed. No E i/DV is included, since the amount of none of the three compounds tracked (FA, CAF, CHA) was higher than the chromatographic detection limit, which indicates that they are mostly found in PV.
Only E CHA/PV , E FA/PV , and E ALL/PV were correctly adjusted to Equation (8). The fitting coefficients of this equation are given in Table 4, where it can be seen that E CHA/PV depends on all of the terms except the cross term (β 12 ), and all of them are significant (p < 0.05) except for the linear term of the CO 2 flow rate (β 2 ). E FA/PV and E ALL/PV also depend on all terms except for the cross term (β 12 ), but only β 0 and the linear and quadratic pressure terms (β 1 , β 11 ) are significant. Figure 5a-c shows the contour plot corresponding to the surface defined by Equation (8) for the E CHA/PV , E FA/PV , and E ALL/PV enrichment ratios. As Figure 5a initially shows, E CHA/PV increases as pressure increases and then decreases at higher pressure values. For a fixed pressure, E CHA/PV behaves similarly with increasing CO 2 flow rate. Then, the maximum is located in a central area, with pressure values between 103 and 150 bar and CO 2 flow rates of 12-53 g/min. E FA/PV and E ALL have a similar behaviour, as In relation to the physical-mathematical treatment of the considered SAF process, it can be indicated that there are close similarities with the well-known supercritical antisolvent precipitation (SAS) processes suitable for the production of fine powders, or even composites including a bioactive substance and a polymer. However, the SAF modelling turns out to be much more complex, because mixtures consisting of a high number of components, such as the natural extracts studied in this work, are considered. Tentatively, most of the developments for the different steps implied in SAS could be adapted to SAF [32]. There are specialized papers about SAS that model both the whole process [67][68][69] as well as specific steps, such as breakup of the feed solution jet in the supercritical CO2 stream [70][71][72], mass transfer between droplets of organic solution and the compressed antisolvent [73][74][75], supersaturation and nucleation rate [76][77][78], growth mechanism [79], or the influence of the variables on precipitation [80]. However, as explained in a previous work [50], detailed modelling of the SAF process is quite a complex task, because of the high number of components in the extracts that usually need to be processed, then becoming rather inextricable from a rigorous theoretical study. In addition, reaction or association processes between the different components can occur. Therefore, and to the best of our knowledge, to date, there are no publications on the detailed theoretical modelling of SAF. In relation to the physical-mathematical treatment of the considered SAF process, it can be indicated that there are close similarities with the well-known supercritical antisolvent precipitation (SAS) processes suitable for the production of fine powders, or even composites including a bioactive substance and a polymer. However, the SAF modelling turns out to be much more complex, because mixtures consisting of a high number of components, such as the natural extracts studied in this work, are considered. Tentatively, most of the developments for the different steps implied in SAS could be adapted to SAF [32]. There are specialized papers about SAS that model both the whole process [67][68][69] as well as specific steps, such as breakup of the feed solution jet in the supercritical CO 2 stream [70][71][72], mass transfer between droplets of organic solution and the compressed antisolvent [73][74][75], supersaturation and nucleation rate [76][77][78], growth mechanism [79], or the influence of the variables on precipitation [80]. However, as explained in a previous work [50], detailed modelling of the SAF process is quite a complex task, because of the high number of components in the extracts that usually need to be processed, then becoming rather inextricable from a rigorous theoretical study. In addition, reaction or association processes between the different components can occur. Therefore, and to the best of our knowledge, to date, there are no publications on the detailed theoretical modelling of SAF.
In any case, it seems clear that the different solubilities shown by the components of the mixture to be separated in the supercritical ethanol-CO 2 mixture play an essential role in the fractionation process. In fact, it can be mentioned here that a very simple semi-empirical model containing the operating parameters (temperature and pressure) and the Hansen solubility parameters has been proved to be appropriate for describing the selectivity in a GAS-fractionation process closely related to SAF and SAS [81,82]. Related to the system considered in this work, it can be pointed out that the feed solutions present a glutinous character, which may be due to components such as resins or mucilage [83,84] not being completely solubilized. The presence of this solid material could induce a very effective process of heterogeneous nucleation [76]. This is how a practically complete precipitation in PV of the widely supersaturated CAF, CHA and FA occurs.
The optimal conditions for optimum yield and enrichment values, indicated above (153 bar, 42 g/min), are located between the extreme values of the ranges studied. This could be explained by the opposing effects that take place both for the increase in CO 2 flow rate and pressure. An increase of CO 2 flow rate implies a higher dragging effect, but also favours better mixing and nucleation. For its part, an increase in pressure would increase solubility, but favours nucleation of solutes.

Results of the Skin Model
The skin model applied was a fully hydrated skin model, and all the calculations for resistances, possible pathways and overall permeability values were performed using the COSMOtherm software package. The predicted values for Equations (9)-(11) for individual compounds are gathered in Table 5.
The compartment with the largest resistance for CAF and FA is the stratum corneum (logR SC = −5.61 and −5.51 respectively), while for CHA, the largest resistance is the stratum spinosum (logR SS = −7.06). This means that while for CAF and FA the stratum corneum is the main barrier in the penetration process, for CHA it manages to penetrate up to the stratum spinosum. Regarding the penetration pathway, both CAF and FA prefer the trans-corneocyte route (R SC, inter >> R SC, trans ), while CHA prefers the transcellular route of the Stratum spinosum, called the interstitial matrix (R SC, inter << R SC, trans ). The overall resistances, R skin , including the shunt pathway, are R skin = R skin = 7.74 × 10 7 , 1.86 × 10 9 and 3.74 × 10 7 for CAF, CHA and FA, respectively. Translating them using Equation (9), the permeability values for each compound were logK p = −5.89 for CAF, logK p = −7.27 for CHA and logK p = −5.57 for FA. These values are corrected by a constant offset (∆logK p = −1.12 cm/s) [46], and the final permeability values are logK p = −7.01, −8.39 and −6.69 for CAF, CHA and FA, respectively. Table 5 also shows the average calculated permeability values of CAF and FA in skin by Zhang et al. [85], and experimental permeability values from Caco−2 cells for CAF, CHA and FA [86,87], as well as the deviation of our corrected values logK p . It can be seen that our logK p corrected values fit better to logK p (calc) of skin by Zhang et al. (deviation = 0.16 and −0.46 for CAF and FA) than to the values of logK p (exp) of Caco−2 cells (deviation = 1.17, 2.79 and 1.71 for CAF, CHA and FA respectively). This may be due to the fact that Caco−2 cell lines are used as a model of intestinal absorption and the estimation of oral bioavailability [88], and therefore, their permeability differs from the permeability of the skin. For a compound to be considered toxic, it must reach a certain layer of the skin. Stratum corneum, consisting mainly of dead cells, is considered a safe layer in which toxicity or irritation does not occur [89]. This makes both CAF and FA safe for cosmetic use. CHA, on the other hand, reaches a deeper layer of the epidermis, but does not reach the dermis, and furthermore, its permeability coefficient is lower than the desquamation rate coefficient (logK D = −9) [90], which also makes it safe for cosmetic use. This desquamation rate coefficient indicates that even if a compound penetrates beyond the stratum corneum, the skin regenerates and prevents the compound from penetrating deeper in the epidermis [89].

Conclusions
The combination of two advanced separation techniques applied to Calendula officinalis flowers was successfully used to obtain extracts of interest in cosmetics. The initial supercritical extraction with CO 2 (SCE) carried out, the yield of which was 8.3% and the subsequent maceration in ethanol applied, the yield of which was 8%, showed results similar to those obtained in works by other authors. Subsequently, the ability of the supercritical antisolvent fractionation (SAF) technique to fractionate the ethanolic extract to obtain an enriched fraction in the tracked compounds (CHA, FA, CAF) was satisfactorily evaluated through a Response Surface Methodology (RSM) based on Central Composite Design (CCD). Chlorogenic acid, ferulic acid and caffeic acid were obtained almost exclusively in the precipitated fraction, thus obtaining a fraction enriched in antioxidant compounds of great interest in the cosmetic industry and even with possible applications in the pharmaceutical or food industry. The CCD design also made it possible to deduce the pressure and CO 2 flow rate conditions to obtain optimal yields and enrichment values were 153 bar and 42 g/min (composite desirability = 0.673).
Regarding the skin permeability model used, it was possible to evaluate the computational permeability of CAF, CHA and FA. In all cases, none of the compounds pass through the epidermis completely, which means that they are safe compounds for topical use on the skin. The final permeability values predicted by the model were logK p = −7.01, −8.39 and −6.69 for CAF, CHA and FA, respectively and were compared with experimental values for both abdominal skin cells and Caco−2 intestinal cells, concluding that the skin model is a good model for predicting permeabilities and behaviours of compounds in the epidermis.
Although the SCE was carried out to favour obtaining the antioxidant phase, the extracts obtained through this technique are of great value due to their composition, and they are widely used in the cosmetic industry. In fact, a cosmetic formulation has been prepared containing extracts of C. officinalis, obtained under the conditions used in this work, that has been approved by the French agency ASNM (L'Agence nationale de sécurité du médicament et des produits de santé), and it is now being tested on human volunteers at MEDES (Institut de Médicine et de Physiologie Spatiales, Toulouse, France). All this makes these two advanced separation techniques (SCE and SAF), together with the experimental design and the skin permeability model used, powerful tools for obtaining, concentrating and evaluating compounds of interest to the cosmetic industry.